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Abstract. We study a model of self propelled particles exhibiting run-and-tumble 
dynamics on lattice. This non-Brownian diffusion is characterised by a random walk 
with a finite; pc;rsistcncc length between changes of direction, and is inspired by the 
motion of bacteria such as E. coli. By defining a class of models with multiple 
species of particle and transmutation between species we can recreate such dynamics. 
These models admit exact analytical results whilst also forming a counterpart to 
previous continuum models of run-and-tumble dynamics. We solve the externally 
driven non-interacting and zero-range versions of the model exactly and utilise a field 
theoretic approach to derive the continuum fluctuating hydrodynamics for more general 
interactions. We make contact with prior approaches to run-and-tumble dynamics off 
lattice and determine the steady state and linear stability for a class of crowding 
interactions, where the jump rate decreases as density increases. In addition to its 
interest from the perspective of nonequilibrium statistical mechanics, this lattice model 
constitutes and efficient tool to simulate a class of interacting run-and-tumble models 
relevant to bacterial motion, so long as certain conditions (that we derive) are met. 
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1. Introduction 

The last decade has seen a growing number of studies of biological systems conducted 
by physicists. Much of this work relies on tools not traditionally found in biology. For 
instance, recent methods of nonlinear optics have made possible observation of biological 
systems on much smaller scales than was previously possible [l]. In other cases biological 
systems have helped shed light on questions of fundamental importance in theoretical 
physics. Studies of bird flocks showed that it may be possible to observe long range order 
in systems with a continuous symmetry if detailed balance is broken (in equilibrium 
this is forbidden by the Mermin- Wagner theorem js]). 

Bacteria provide one example of a biological system which is of great interest to 
the study of nonequilibrium physics. Some species of bacteria, such as Escherichia coli, 
swim by means of a series of relatively straight runs interspersed with short periods of 
chaotic motion, known as tumbles, during which their orientation changes at random and 
they experience little net displacement |4] . This can be seen as a type of non-Brownian 
diffusion where the steady state probability distribution will not be Boltzmann. 

Following experiments to determine their behaviour by Berg and others [5 11 



since the 1970s, much is known about the dynamics and behaviour of individual 
bacteria. Less, however, is known about their collective behaviour and it is here 
that statistical mechanics can play a useful role. Schnitzer et al. analysed various 
strategies bacteria may employ, e.g. changing their speed or tumble rate, and studied 



the differences these make to the steady state distribution 12 . These analyses, which 



were based on off-lattice descriptions of the particle dynamics, were more recently 



extended to allow for various types of interaction between bacteria 13 , external 



force fields 14 , hydrodynamic interactions 15 , and coupling to logistic population 
growth 16 . Intriguingly it was shown that a dependence of motility parameters such 
as swim-speed on density was alone enough to cause phase separation 13 , in contrast 
to detailed-balance systems where many-body effects on particle mobility cancel exactly 
in Boltzmann steady states. 

In this work we construct a lattice model whose non-interacting continuum limit 
recovers the off lattice equations previously derived for noninteracting run-and-tumble 



particles 12 , 14 . We will also address interactions in the form of density dependent 



motility parameters, and thereby connect also with the off lattice approach of 13 
One motivation for this approach is as follows. Although real bacterial systems are 
off lattice, both the microscopic run-and-tumble equations, and the diffusive continuum 
equations found by coarse graining these, are difficult to simulate efficiently, particularly 
once interactions are included. On lattice the local density of particles is easy to 
determine and so density-dependent interactions are easy to include at relatively low 
cost computationally. Lattice simulations are also easier to extend to higher dimensions 
as we will consider towards the end of this paper. For all these reasons, creation of 
a robust and accurate lattice representation of bacterial motility is a worthwhile goal, 
even from a purely phenomenological standpoint, whereby the purpose of a model is to 
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provide a fairly direct explanation for results seen in experiment. 

A second motivation is more fundamental. Models of non-equilibrium statistical 
mechanics can be broadly split into two categories: one phenomenological as just 
described, the other comprising simple models which allow the study of basic concepts 
and facilitate a more detailed understanding of the nature of non-equilibrium physics. In 
the latter category we can think of lattice transport models such as the exclusion 17 or 



zero-range 18 processes, for which some exact analytical results can be found, as well as 



methods to characterise fluctuations and large deviations in non-equilibrium states 19 



Our model falls squarely into this category of simple theoretical models, and indeed 
extends some of these examples; we investigate both a zero-range process and a partially 
excluding version of our run-and-tumble system. As a microscopic lattice model we can, 
under certain conditions, calculate exact steady states and understand precisely how 
changes in the underlying dynamics affect the probability distributions. More generally 
we can always write an exact master equation for probability of a configuration and 
utilise a variety of field theoretic representations to derive the large scale fluctuating 
hydrodynamics, to attempt to map the system to a free energy and to determine the 
steady state behaviour and dynamic stability. 

As already explained, in addition to being an interesting toy model that we can 
use to understand driven non-equilibrium physics, our model also forms a direct lattice 
counterpart to some well established phenomenological continuum models of bacterial 
dynamics. It is unusual for a model to allow exact computations while credibly 
describing the real behaviour of a physical system, and this is one of the reasons why 
we consider this model to be of interest. 

In section [2] we introduce our model of run-and-tumble bacteria on a lattice. 
In section [3] we consider those cases in which we can find exact solutions to the 
steady state behaviour. After first examining in section 3.1 the non-interacting case, 
which for inhomogeneous and anisotropic jump and tumble rates can still admit non- 
trivial solutions, in section 3^ we include a zero-range interaction and derive sufficient 
conditions on the rates for a factorised steady state to exist. 

In section |4] we consider more general interactions, where the jump rates can 
depend on the full configuration of the system, and derive a continuum fiuctuating 
hydrodynamics using a field theoretic approach. We first set up the theoretical apparatus 
and re-derive the results for non-interacting systems within this framework; we then 
calculate the field equations for the more generic interacting system. In section [5] we 
consider one specific type of interaction and investigate the effect of crowding. We derive 
a free-energy-like functional from which we can determine the steady state probability 
distribution and compare the effects of different methods to coarse grain the local density 
in the interaction terms. Finally we extend our results to higher dimensions more 



applicable to real systems in section 5.3 
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Figure 1. Presentation of the model. Filled circles represent right moving particles 
while unfilled circles denote left moving particles. Some of the possible transitions are 
illustrated on the figure. 



2. The Model 



To model the finite persistence length in run-and-tumble dynamics on lattice particles 
jump repeatedly in the same direction, u, with rate d{\i) and change direction - tumble - 
with rate a(u). In Id, this means particles can be either right-going or left-going. Right- 
going particles jump with rate df from site i to site i + 1 and tumble with rate af . After 
a tumble, they become left-going particles with probability 1/2. The corresponding rates 
for left-going particles are called d'^ and a^" , see figure [l] Throughout this paper we 
assume tumbles to be instantaneous. In reality the duration of tumbles is typically of 
the order of one tenth of the duration of runs |4j. There may be situations, however, 
where time spent tumbling may increase, and where the finite duration of a tumble may 
have an effect on the dynamics and steady state of the system. We leave a consideration 
of this case to future work. 



3. Exactly Solvable Cases 

The jump and tumble rates defined in section [2] may, in general, depend on the 
occupation numbers of any lattices sites. Through this we may introduce interactions 
to our system and consider the general case of interacting bacteria. Before we deal 
with such complex cases in section |4[ however, let us look at the interesting limiting 
cases that allow for exact computation of the steady state. In particular we start by 
considering first the non-interacting case and investigating the effect of position and 



direction dependence on the rates in section 3.1, before then turning to zero-range 



processes in section 3.2 



3.1. N on- Interacting Particles 

For non-interacting particles, one can always handle the single particle case first and 
compute the average occupancy pf of left (— ) or right moving particles on site i. 
As shown below, the steady-state distribution of n particles is then given by a product 
on each site of a multinomial distribution of parameters pf . 

Let us start with the single particle process and call P{k., +) and P(A;, — ) the 
probability to find a bacterium at site k going to the right and to the left, respectively. 
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The master equation reads 



d,P{k, +) = dUP{k - 1, +) - dtP{K +) + -) - +) (1) 

d.Pik, -) = d,^,Pik + 1, -) - d-,Pik, -) - -fPik, -) + -fPik, +). (2) 

3.1.1. Exact Solution of the Steady State For the case of non-interacting particles, 
where all interesting and non-trivial effects come from heterogeneities or anisotropics in 
the jump and tumble rates, we can calculate the steady state probability distributions 
exactly. We start from the continuity equation 

t[P{t, +) + P{t, -)] = Ji_i,i - Ji,,+i, (3) 

where Ji,i+i is the net probability flux between sites i and i + 1, that is 

J,,i+i = di P{i, +) - rfr^i Pit + 1, -). (4) 

As we are only after the steady state distribution we take the time derivative on the 
left hand side of equation ^ to be zero, so all currents Ji,i+i are constant and equal to 
some J fixed by the boundary conditions. 

The master equation for the density of left-moving particles at site i reads 

= dr^, P{z + 1, -) - dr P(z, -) + ^P(^, +) - ^P{t, -). (5) 

Using Q to eliminate the first term on the right hand side of equation ([s]), which is the 
only term to depend on site i + 1, we can establish a relationship between the number 
of left and right moving particles on site i at steady state: 

Pih +) (dt + ^) - P{h -) (d. +^)=J- (6) 

Equation (|6|, along with equation (|4]), leads to the recursion relation for right moving 
particles 

Pi^ + 1, +) = iL 2^m + «m p(,^ _ J f.i 
which can be solved to yield 



i-l 



-jy '''''' . . r) (8) 



2d- +a. 

The probability to find a particle at any position and in either state can then 
be calculated by noting that the total distribution must be normalised, i.e. 
[P{i, +) + P{i, —)] = 1 and that J is imposed by the boundary conditions. For 
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example, closed boundaries require that J = 0, while for periodic boundaries we have 
the additional constraint that P(L + 1, ±) = P(l, ±). Note that the probability densities 
for left and right moving particles do not have to be the same, and, in general, will not 
be. 

The probability of a given configuration P(|n^, }), is given by 

where N is the total number of particles in the system. If we call the average number of 
right or left going particles on a site pf = N P(i, ±), then in the limit where N ^ oo and 
L — )■ cxD, so that P(z, ±) — )■ 0, but pf remains finite, the probability of a configuration 
is given by 

j ■ ■ 

3.1.2. Continuous limit Since we ultimately want to compare the run-and-tumble on 
lattice with its off lattice counterpart, let us first take the continuum limit of the master 
equation. Explicitly introducing the lattice spacing a and defining Xk = ka, the master 
equation ^ reads 

dtP{xk, +) = - a)P{xk - a, +) - d^{xk)P{xk, +) 

+ ^Pfe.-)-^^P(...+) (12) 
dtP{xk, -) = d'{xk + a)P{xk + a, -) - d'{xk)P{xk, -) 

- ^^P{xk, -) + ^^P{xk, +). (13) 

We are interested in cases where the typical run length is much longer that the lattice 
spacing so that rf^ ^ a^. Furthermore, when a — ?■ while f^(a;) {x)a remains 
finite, one gets, at leading order, 

dtP{x, +)= - V[t;+(x)P(x, +)] + ^^P(x, -) - ^^P(x, +) + 0{a) 
d,P{x, -) = V[v-{x)P{x, -)] - ^P(x, -) + ^^P{x, +) + 0{a), 

(14) 

which is exactly the master equation for run-and-tumble bacteria considered previously 



off lattice 12 , 13 . Following the same path as there would lead to a Langevin equation 
for the density of a large but finite number of bacteria 

p{x,t) = -V[pV -DVp+^/2D^r]], (15) 

where 
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with a = (a+ + a')/2, v = (t>+ + v~)/2. It was shown that ( |15| ) captures the steady 
state of the off lattice model exactly |13j and in section |4] we will show that it also 
describes the large scale behaviour of run-and-tumble bacteria on lattice. The condition 



for (15) to derive from an effective free energy is that there exists an excess free energy 



functional J-'ex [p] such that 



D 6p{x) 
which can be solved to give 

-L 



V^-^, (17) 



J^ex[p\ = J dx <P{X) log J + o / 



2 /n \V+ V 



as long as, for periodic boundary conditions, dx{a — a^v ) = 0. The total free 
energy is then given by 

J^[p] = [ dxp{\ogp - 1) + J^M. (19) 



Note that it is not always possible to write a free energy of this form for non-interacting 
particles in higher dimensions, nor, in general, for interacting systems. As we can see, 
there is no gradient term in this expression. This is due to the fact that when deriving 



the continuum limit, terms of order a and higher are neglected. If equation (15) leads 
to large gradients in the density, these higher order terms should be included and may 
alter the result; these terms control the surface tension, for example. These higher order 



terms could also violate the condition (17). 



We will now illustrate the steady state predictions of equations ([s]) ^ and (19) for 
a few simple cases and show the difference between the results predicted by the lattice 
model and by the continuum approximation. 

3.1.3. Examples Where not stated otherwise the simulations we present here use 
reflecting boundary conditions; if a particle tries to jump off one end of the lattice 
it is, instead, kept where it is but turned around. All simulations are performed with 
continuous-time Monte Carlo algorithms. 

Position Dependent Rates with Closed Boundaries First, we consider the case of a 
position dependent, but isotropic, jump rate and a constant tumbling rate. As a simple 
example we use a top-hat function for jump rate such that = 1+10 ^(z-150) e(350-i), 
where 6{x) is the Heaviside step function. Both the continuum and lattice theory predict 
that the average occupancy should be inversely proportional to the velocity, 

1 _ 1 

p{x) oc — ; p, = pf + Pi oc -. (20) 



The results of the simulations and both predictions are shown in figure |2] (main). 

In contrast to the jump rate, simply making the tumble rate depend on position 
but maintaining isotropy has no effect on the predicted distribution. Note that the free 



energy in equations (18) and (19) has no dependency on a for isotopic rates. Using 
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Figure 2. Main: Steady state probability distribution for constant tumble rates, 
of = 1 and isotropic jump rates df = 1 + 10 9{i — 150) 0(350 — i). Data averaged from 
the positions of 400, 000 particles. Inset: Steady state probability distribution for 
constant jump rates, df = 10, and tumble rates af = 6{i — 150) 6'(350 — i). Data from 
100, 000 particles. In both figures simulation data are shown in red and the theory 
prediction in blue. Both simulations performed on a lattice of 500 sites and recorded 
at t = 5000. 



the same form as for the position dependent jump rate in our simulations this can be 
verified numerically, see figure [2] (inset). 



Direction Dependent Rates with Closed Boundaries In many physical situations, 
however, bacteria do not move unbiasedly but are affected by their external conditions. 
This may be due, for example, to sedimentation due to gravity, where there is an 
asymmetry in jump rates between left and right (or up and down) moving particles. 
Another case of interest may be anisotropic tumble rates. Bacteria undergoing 
chemotaxis often vary their tumble rate dependent on whether they are travelling up or 
down a chemical gradient. Though a simple asymmetry in tumble rate does not fully 
capture this behaviour, we do see particles preferentially move in the direction of a lower 
tumble rate, as would be expected. 

These two cases show qualitatively the same behaviour, with an exponential decay 
in the unfavoured direction. From equations (|8]) and (|9]), the probabilities for left or 
right going particles on lattice are 

'd+ {2d- 



PH,±] 



oc exp 



i log 



a 



exp [i Aiatt] 



d- {2d+ + a+ 

while in the continuum case the probability density is given by 



p(x) oc exp 



a V 



a ' V 



X- 



exp [x A 



cont I 



(21) 



(22) 



2v+v- 

For expositional simplicity we consider homogeneous rates here. To examine the 
difference between these two cases, consider, for example, the case of sedimentation, 
where o?^ = (io(l ± e) and = olq. The decay constant in the continuum limit is then 
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Figure 3. Steady state probability distribution for constant tumble rates, af = 1 
and jump rates df = 10^ 1. Simulation data are shown in red, the lattice prediction 
in blue and the continuum prediction in pink. Data collected from 10^ particles at 
t — 2000 on a lattice of 500 sites. 



and the lattice decay constant is given by 

We can see then that the two decay lengths will be equal if the jump rate is much larger 
than the tumble rate, i.e. for large average run lengths. Both decay constants tend to 
zero as the asymmetry disappears, e — )■ 0, but the ratio Aiatt/^cont remains finite. 

In our simulations we use df = 10^1 and ao = 1, so the drift velocity, the external 
bias, is much less than the self-propelled speed, i.e. \d'l — d^\ <^ c?^. The continuum 
theory predicts the distribution to be p{x) oc exp(— x/99), while the lattice theory 
predicts P{i) oc (207/209)*. Both predictions are shown, along with the simulation data 
in figure [3] The ratio of decay constants is 

1 - -^j"^ + o(^X^ 0.95. (25) 



Acont do(l-e^) \do 



Note that in equations ([23j) and (24) one constant multiplies lattice position, i, and 
the other the continuum position, x, hence the factor of a. We see that the difference 
between the lattice and continuum results vanishes in the infinite run length limit, 
do/ao — 7- oo, unless e = 1, in which case both decay constants diverge. 

Periodic Boundary Conditions We can also calculate the expected probability 
distribution for periodic boundary conditions. In this case our calculation on lattice 
is slightly more complicated as we do not know the current a priori, but must determine 
it through the conditions P{L + 1, ±) = P(l, ±) and J2i +) + Pih -)] = 1- 

We can write P{i, +) = Ci{i)[P{l, +) — JC2(i)], where Ci{i) and C2{i) can be read 
from equation ([8|. Then from the periodicity of the system we can write the current as 
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Figure 4. The probability distributions at steady state for rates given by 
equations ( 28 ) . The distribution for left moving particles is shown in blue and for 



right moving particles in red. The points show data from stochastic simulations and 
the solid lines the theoretical prediction. Data from 2, 000, 000 particles at t = 200 on 
a lattice of 200 sites. 



We can then use the normahsation of the distribution to determine -P(l, +) as 



T- (27) 

As an example see figure |4] where we consider the case where the jump and tumble 
rates are 

dt = 10 d7 - ^ 



exp(-(a;-100)2/5000) 

at = 1 a7 = 1. 



(28) 



We omit the exact forms for the probability distributions as these do not reduce to a 
compact form and are not enlightening in themselves. 



3.2. A Zero-Range Interaction 

Though there is no generic solution for the steady state of interacting run-and-tumble 
particles, there are limiting cases that can be solved exactly, most simply a zero range 
interaction. We define the number of left and right going particles on a site as, 
respectively, n~ and nf. The total occupancy of a site is then = nf + n~ . The 
simplest interaction we can add is a zero-range interaction, where the rate for a particle 
to jump from site i with occupancy (nf^n^) to site i ± 1 is defined as df{nf,n~) and 
is a function of the number of particles at the departure site but not dependent on 
the number of particles at the arrival site. With this addition we can now see more 
complex behaviour and non-trivial steady states, even for homogeneous and isotropic 
jump and tumble rates, but can still, under certain conditions presented below, calculate 
the stationary probability distribution exactly. 

As is standard for zero-range processes (ZRPs) [18j, we begin by assuming there 
exists a factorised form for the steady state probability distribution of the form 

L 

P{{nt,nr})cxl[g,{nl,n-). (29) 
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This ansatz can then be substituted into the master equation for P (jn^j-n,"}), 

L 

dtP = + 1) dtint + l,n-,) Pint + 1, n^^, - 1) - n+ n^) P(n+, n^^,) 



k=l 



+ ^(nt + l)P{n, - l,n+ + 1) - ^n+PK,n+). (30) 

Then, for periodic boundary conditions, one way we may choose to solve this equation 
is to separately balance the fluxes for right moving particles entering and exiting each 
site, left moving particles entering and exiting each site and particles tumbling between 
species on each site. We then arrive at three sufficient conditions on the allowed rates 
for such a factorised form to exist: 

gi{n^, n^ ~ 1) = en" d^{n'^, n~) gi{n^ , n^) (31) 
gii^n'^ — l,n^) = c n'^ d^in^ , n^) Qiin^ , n^) (32) 
gi{n^,n^)n'' a~{n'^ ,n~) = gi{n^ + l,n^ - 1) (n+ + 1) aKn"^ + l,n^ - 1), (33) 

in which c and c' are arbitrary constants. The first two of these conditions are the same 



as Evans and Hanney found for their two species model without transmutation 18 , 20 



while the third is due to balancing the tumbling. Note that, in principle, there may be 
other ways in which we can balance these terms which could arrive at different conditions 
on the rates. 



Putting the three conditions (31)-(33) together and eliminating the factors 
gi{n~^,n~) yields two constraints on our choice of rates: 

d; (n+,n-) dl (n+,n" - l) = d; (n+ - 1,71") rf+ (34) 
n^ al {n^,n~) _ c {n^ + l)di {n~^ — l,n~ + 1) ^^^^ 

{n~ + l)a~ (n+ — 1, n~ + 1) d n+ df (n+, n~) 

One natural, but again not necessary, way to fulfil these conditions is to take 

df{n^, n~) = uf{n~^) ujii^n'^ + n") 

d~ {n^ ,n~) = (n^) uji{n'^ + n^) (36) 

for the jump rates. That is, for the rate at which a left or right oriented particle moves 
to be a product of a function of the number of particles oriented in that direction, and a 
function of the total number of particles on a site. Both functions can vary from site to 
site; the first can also depend on the particle species, but the second must be the same 
for both. 

A sufficient conditions on the tumble rates is then to take 

at{n'^,n'') = cut{n^)Ai{n'^ + n^) (37) 
a~ (ji^ ,n~) = cuj {rT ) Ai {n^ + n~). (38) 
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Figure 5. Left: The effective free energy density for tfie zero-range interaction witfi 
jump rate given by equation (41 1 for af = 1, {rii) = 12, — 20 and Vq = 10. Right: 
A typical snapshot of the system during its relaxation towards equilibrium for the same 
parameters on a lattice of 200 sites and 2400 particles 'di t — 1000. 



The functions u'l{n^) and u~{n~) are the same as in (36 ); Ai{r) is a new, unconstrained, 
function that appears in both rates. 

Up to a constant that can be subsumed into the normahsation, the marginals are 
given as 

n+ ^ n~ ^ " 1 1 "~ 1 

,.in\n-) = n-57(^ n^fM = n-^ n^. (39) 

where 7 = c/c' and = + n~ . We can then re-write the probabihty of a given 
configuration as 

Pii^^^i}) = ^tl9^{nt,n-) = leSiiK^^(""'"n), (40) 

where Z is a normahsation, and define an effective single site free energy fi{nf,n~) = 
— \og{gi{nf,n~)). Note that this is independent of af{n~^,n~); the way in which we 
chose to solve the master equation does not lead to factorised steady states for which 
the distribution can depend on the tumbling rates. As we saw that asymmetric tumble 
rates could affect the equilibrium distribution for the non-interacting case, we might 
suppose there are other solutions for the zero-range process which admit distributions 
dependent on the tumble rates. Whether or not these allow for factorised steady states 
remains to be determined. 

To foreshadow the finite range interaction we will examine in section |4| and to 
mimic the situation where an increase in density decreases the particles motility (as, for 
example, they get in each other's way) we now consider the following particular form of 
the steady state for this two-species ZRP for jump rates 

,±/ + J Vo[l - {n+ + n-) /rim] if n+ + < ra^ 

dfin^.n ) = \ I r 4. _ (41) 

y vo/nm itrT + n > Um 

and tumble rate af = a. That is, the tumble rate is constant per particle and the jump 
rate decreases linearly as density increases until reaching a constant rate of Vo/rim at 
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Figure 6. A time averaged snapshot of the steady state of zero-range process with 
parameters Um — 20, {n) = 12 = 9, = 11, af = 1. Data averaged from 10,000 
snapshots between t = 40, 000 and t = 50, 000 



n'^ + n~ = rim — 1. In this case the effective free energy is double welled and the system 
separates into isolated sites of high and low density. The relative numbers of high and 
low density sites to which the system first separates are initial condition dependent. The 
system then relaxes via a series of evaporations and condensations towards a fixed steady 
state. This is reminiscent of what happens for single-species zero-range processes with 
similar jump rates. Where the jump rates remain finite even for very large occupancies 



the normal condensation is arrested 21 . The free energy is shown in figure [5j along with 
a typical snapshot of the system. 

To this behaviour we can then add a drift term to simulate sedimentation by biasing 
the jump rates in one direction and applying closed boundary conditions. We see all 
the high density sites collect at one end of the lattice and the low density sites at the 
other, see figure [61 



4. Fluctuating Hydrodynamics for Interacting Bacteria 

We now turn to the more general case of interacting bacteria with interactions that 
are not limited to being on-site. Specifically, we allow the jump and tumble rates to 
depend on the occupation numbers of each lattice site so that 

df = df{nf)- af=af{nf), (42) 

where nf is a coarse-grained occupancy that depends linearly but non-locally on the 
occupancies of the whole lattice, 

^? = E^-.^r (43) 
j 

In general the coarse graining kernel could also be a function of lattice position 
though here we do not consider that situation, where the manner in which the density 
is felt by the particles varies with position. We aim at describing large scale behaviour, 
i.e. on a colony size, and so to derive a fluctuating hydrodynamics. In addition, this will 
allow us to compare again the phenomenology on and off lattice and to look for cases in 
which there is a "free-energy" like description, and for which we can thus characterise 
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the steady state. We follow a field theoretic approach to derive a continuum Langevin 
equation for the system, from which we can deduce the appropriate Fokker-Planck 
equation and the steady-state distribution. 



4.1. Field Theory for Non-interacting Particles 



To illustrate the technique we shall use to construct the fluctuating hydrodynamics 
for the full interacting case, let us first handle the non-interacting case and re-derive 



equation (15). 



Field theoretic representations of lattice gases using bosonic coherent states were 

The case where each site 



established in the 1970s following Doi and Peliti 22 , 23 



is limited to a single particle can be handled in some cases in this formalism 24 



while more general finite occupancies could be handled using spin coherent states 25 



Alternatively, probabilistic approaches from mathematical physics have also been 



used 26 -28 . Here we use an alternative derivation, based on an approach a la Jansen 



and De Dominicis [29 30 transposed in the context of the master equation. This is very 
similar to the generating function approach used by Biroli and Lefevre 31 . Beginning 



with a process discrete in both time and space one writes the probability of a trajectory 



as 



N 



P[{n+(t,),nr(t^.)}] = {\{\{Knt{t,^i)-nt{t,) - Jtih)) 

i=i j=i 

X 5(nr(t,+i) - n-(t,) - J-(t,)))_^, (44) 

where nf{tj) is the number of right (+) of left (— ) moving particles at site i at time tj, 
the Ji{tj)^ are the changes in the number of each type of particle at each site at each 
time step and the bold faced J denotes the average is over all J's. Re-writing the Dirac 
delta functions using imaginary Fourier representations this can be written as 

P[{nt{t,),n;{tj)}] = / ]J]Jdht{t,)dhr{t,)(exp(nt{t,){nt{t,+i) - nt{t,) - J+(t,)) 

i=i j=i 

+ niit,)inr{tj^,) - nr{t,) - ^r(i.))))_j, (45) 

where it should be noted that the conjugate fields hf are imaginary. The average over 
the J's can then be calculated explicitly from the dynamics. Specifically, a right moving 
particle can jump from site i to site i+l at time step j with probability nf{tj)dfdt, where 
dt is the duration of a time step. The corresponding values of the J's are Ji(tj) = — 1 
and Ji^iitj) = 1. Calculating all other moves and the probability that nothing happens, 
which corresponds to all jf = 0, we can write 

(e-^<^^)Jt(ti)-K(t.W(t.)'j_^ = l + ni{t,)+dt (^e"»+fe)-"iife) - l) dt 



+ (e"^(*^^-"'"(*^^ - 1) dt 



dt 
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+ ^^r(^j) (e"^"(*^^"">^(*^^ - l) dt. (46) 

As this is of the form 1 + A;dt we can approximate it as exp(A;dt) and write the probabihty 
for the trajectory as 

^ . L N V r- L N 

P[{nt{t,),ni{t,)}]= [l[l[dht{t,)dni{tMexp 



n. 



a,, 
+ —n 
2 



a, 
+ —n 
2 



We can then take a continuous time hmit and make the substitutions 

N ^T=Ndt N 



nf{t j+i) - nf{tj) ^ n^^dt; XI ^ 







dt; l[dnf{t,)-^Vnf. 



The probabihty of a trajectory can then be written 



P[{ntit),ntit)}]= [HV 
where the action S is given by 



(47) 



(4J 



(49) 



S 



dt J2 [nffif + fiih- + ntdt (e^^-"ii - l) + nr^ic/'+i (e^^'+i""'" - 1 



1 



(50) 



Note that generic changes of variables in (49) will result in Jacobians. If these do not 
depend on the fields, nf and hf^ they can be subsumed into the normalisation of the 
path integral but they must be handled with care otherwise. 

We further simplify by considering symmetric, constant, rates d^ = d~ = d and 
af = = a; the more general case causes little conceptual difficulty but is considerably 
more cumbersome as an illustration. Let us then introduce the new variables 

p. = nf + n-; Ji = d{nf-nr)- p. = i(n+ + nr); J. = l(fi+ - n^) (51) 
The action can then be written as 
S= - 







dtY^ \p^P^ + ^Jji + -pi ('e-Cft+l-P' + ^^+l-'^^) + ^P^+l-p^-a+l-J^) _ 2 



Ji 



_j_ _1 I ^--{Pi + I—Pi + Ji+I — Ji) Qpi+l—pi — {J-i 



i + \ — Ji)^ 
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+ f (A+i - P. - ^^^±^) (e''--P-(^--^') - l) 
4 V / 4(i V 



(52) 



The continuum limit can be taken by explicitly introducing the lattice spacing a and 
making the substitutions 



i>t=lja 

Pi — )■ ap{x); Pi — )■ p(x); d — )■ t'a^"'^; / dxa~ 

Jo 



1. 



Ji -> J(x); J, -> J(a;); aV + (53) 

where Vj is the discrete gradient, e.g. ViPi = Pi+i — Pi, and i is the system length. 
After Taylor expanding the action in powers of the lattice spacing and taking a diffusive 
rescaling of time and space, see Appendix A we find that the fluctuating hydrodynamic 
action is given by 

dt j dx [pp - vpVJ - JVp + apP + j , (54) 

which is invariant under a further diffusive rescaling of space and time. 

Going back to the definition of the probability (49), we can then work backwards 
to recover a continuity equation for p from our action [25j. Starting from 

P[{p(a;,t), J(x,t)}] = i ! V[pJ]e~'^^^P^'^P^'\ (55) 



one can remove the quadratic term by introducing a new field ?7(x, t) via a Hubbard- 
Stratonovich transformation so that 

P[{p(a;,t), J(x,t)}] = ij D[p, J,r7]e-^»[^'^'^'^'''l, (56) 
where the new action now reads (after some integration by parts) 

So = - dt dx (pp + vVpJ + pVJ + ^2^pT]J +- 7^r?^).(57) 

Jo Jo ^ V 2 / 

The integral over p and J then leads to 

P[{p{x,t),J{x,t)}]oc y"l?M5(p + VJ)5(^J+ V^^ + Vt;p) e-5/d-dtr,2^ ^^^^ 
where the delta functions impose the two dynamic field equations 

p = -yj- J = -DVp+ j2Dpr]; D = —. (59) 

a 



Given its weight in (|58j), rj{x,t) is a Gaussian white noise: 

(r/(x, t)r]{x', t')) = 6{x - x')6{t - t'). (60) 

This is consistent with the calculation off-lattice for non-interacting, homogeneous and 
isotropic systems and vahdates the results obtained previously, see equation (15). 
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4- 1.1. Fluctuating Hydrodynamics and Large Deviation Functions Before going any 
further, let us make a brief detour to consider the connection with the standard 
fluctuating hydrodynamics approach considered in the mathematics literature [32]. Let 
us first note that from the definition of the continuum limit, one has 

N = y^Pi= [ dxp(x) (61) 
i Jo 

The integral of the density field is thus an extensive variable. On the other hand, the 
density field considered by mathematicians is often defined by 

P(^) = ]j'^Pi^(^ - (62) 

i 

and satisfies the normalisation condition 

N 

dxp{x) = — (63) 

To make the connection between the two approaches, it is thus natural to rescale our 
density term to make the extensivity apparent :p — )• ip. To ensure that the conservation 
equation still has the form p = — VJ, one must also rescale the current field J — )■ £J. 
Before introducing the f]{x,t) field, the action thus reads 

So = -i dt dx (pp - vpVJ - JVp + apJ^ + ^) . (64) 

One can again introduce the noise field and integrate over the conjugates fields p and J 

to get 

P[{p{x,t),J{x,t)}]oc jv[7]\5{p + VJ)5{^J+^2^p7] + Vvp^e-y^''^''^\ (65) 

Interestingly, the fields are now all intensive, and the smallness of the noise does not 
come from a ^/p versus p noise prefactor, but from its explicit variance, read in the 
Gaussian weight: 

{7^{x,t)r^{x' X)) = l5{x-x')5{t-t'). (66) 

This is the usual fluctuating hydrodynamics, as considered, for instance, in In the 
large size limit, the first order correction to the deterministic equation in a l/£ expansion 
is given by the addition of the noise term ^2Dp. This noise is typically of order 1 / 
i.e. trajectories of probability of order 1 have ir]^ of order 1. Large deviations correspond 
to trajectories where the noise can be of order one. They yield probabilities of order 
0(exp(— £)), and are described by the fluctuating hydrodynamics constructed here. 

4.2. Interacting Particles 

Consider now the case of interacting particles where the jump and tumble rates depend 
on the occupation numbers of each lattice site. Our velocity is then modified to 

v^{x) ^ v^{p^{x),x), (67) 
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and the tumble rate to 

a^{x) a^{p^{x),x), (68) 
with p^(x) given by an integral over the density 

p±(x)= / K^{x-y)p{y)dy. (69) 



Following the same path as that followed in section 4.1 for the non-interacting case, 
one gets for the action 

v+{p+)v-{p~) , a+{p+)v'{p~) + a-{p-)v+{p+) ^ 



S = - dtdx 



pp pVJ - JVp + pJ 

v{p+,p ) v{p^,p ) 



+ ———JJ + i pJ 

v{p+,p v{p+,p 



(70) 



where v = {v~^ + v~)/2 and a = {a^ + a~)/2. The factor of i in the final 
term implies that, for the diffusive scaling to hold, at a scale i, the asymmetry 
a~^{p^)v~ {p~) — a^{p^)v^{p'^) must be of order This is reminiscent of the ASEP, 
where if the bias is much smaller than l/y/i the diffusive scaling holds (Edwards- 
Wilkinson universality class), as in the symmetric exclusion process, but for larger 



asymmetries the dynamic exponent z is the same as Kardar-Parisi-Zhang scaling 33 
Integrating over p and J now yields the set of field equations 



p=-VJ; j=-DVp-Vp+\l^ J ^-pv, (71) 



with 



D = -—- V = i— ^^--V^^ (72) 

a 2a; a f 

This formalism can now be used to analyse the effect of interactions on the large 

scale behaviour of a system of run-and-tumble particles. 

5. Crowding Interactions 

Having set up and utilised a field theory apparatus to derive a fluctuating hydrodynamics 
for a general linear dependence on density in the jump rates, let us turn now to a specific 
class of interactions. In particular we consider a crowding interaction, where the velocity 
of the bacteria decreases with increasing density, in which case we expect to see our 
system separate into regions of high and low density as particles become trapped in 
regions of high density jl3j . 

In general we expect to see qualitatively similar behaviour for any choice of f (p) 
which decrease sufficiently quickly towards a finite non-zero velocity at high densities. 
In the following we use 

^ (1 - (p-^)/pJ ifp-^<P. (,3) 

[ ^^o/Pm if P > Pm 

as we did for the exactly solvable zero-range process. 
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5.1. Zero-Range Interactions 



This approach can describe many types of interaction, in particular let us now consider 



the zero range interaction we met in section 3^ in the context of the fluctuating 
hydrodynamics we developed in section |4j This provides us with a benchmark to check 
that our fluctuating hydrodynamics is consistent with the exact results we obtained 
previously. 

For the zero-range process, where the velocity depends only on the occupation at 
the departure site, our kernel is given by Kf = 5iQ and in continuum p^{x) = p{x). 
For simplicity, and to compare with our previous results we shall take a^(p) = a. This 
simplifies equation (71) considerably and, indeed, guarantees that a^v' — a~v^ = 0. 



From section 4.2 we know that the fluctuating hydrodynamics describing the run-and- 



tumble bacteria with velocity given in equation ( 73 ) are given by 



P 



D 



-VJ; 



a 



1 - 



V 



^0 1- 



1- ^ 



a 



The corresponding Fokker-Planck equation is given by 



V 



dx 



5 



5p{x 



a 



1 



_P_ 

Pm 



9. 



P 



p - Dd-^p - Dpidx 



' 6p{x) 



(74) 



V. (75) 



Note that a term V^[-D(p)] could be present, but vanishes for symmetry reasons 



(See 13,34]). Looking for a free energy P oc exp[— J^[p]] one gets 



- V 



6^ 
5p 



-V 



logp + log 1 



P 



whose solution is 
I'[p] = / da;p(logp 



1) - {pm - p) 



log 1 



pr 



for p< pr, 



For p > Pm the free energy density f{p{x)) is given by 



f{p{x)) = p( log 



_P_ 

Pm 



(76) 



(77) 



(75 



which corresponds precisely to the free energy calculated exactly in section 3.2 for the 
total occupancy and in the continuum limit. An example of this free energy for one 
choice of parameters is shown in figure [s] (left). 



5.2. Finite Range Interactions 

We saw in section |5.1| that our fluctuating hydrodynamics admit a free energy and 
accurately reproduced the exact results for the zero-range process. Let us now extend 
that analysis to a system with finite range interactions. As before we now take the 
coarse grained density p"^ to be given by 

p±(x)= [ K^{x-y)p{y)dy. (79) 
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25 50 75 100 125 150 175 200 

Lattice Position 



Figure 7. Snapshot of typical density profiles for an average run length of 100 sites 
(red) and 10 sites (blue) for the isotropic, Gaussian kernel. The black lines show 
the predicted average high and low densities. Data recorded at i = 1000 using 5000 
particles with = 50, k = 2 and a = 1. 



For smooth profiles, we hope that the differences between and p are small so we can 



treat the free energy in equation (77) as a mean field theory. This way we can still use the 



free energy we derived in section 5.1 to predict the coexistence densities and instability 



to spinodal decomposition. The finite range nature of the interactions will introduce 



correlations between sites which we hope will manifest only via surface tension 13 . We 
hope that this surface tension will only effect a clustering of the high and low density 
sites without affecting the coexistence densities. If the mean field theory captures the 
picture correctly, spinodal decomposition occurs whenever the second derivative of the 
free energy density is negative, i.e. for 

- + Vlogt;[p] < ^ — < 0. (80) 

P P Pm- P 

Thus, whenever pm > p > Prnf^, the system should be unstable with respect to spinodal 

decomposition. 

5.2.1. Isotropic Kernels If we use an isotropic kernel to calculate p^ in our jump rates 
we do indeed recover results consistent with the zero-range free energy. That is, if we 
take K^{x — y) = K~{x — y) our simulations match the free energy predictions. In 
particular, we have worked with a Gaussian kernel, where K^{x) = exp{—x'^/k)/Z, 
with Z a normalisation and k a parameter to control the range of the interaction. The 
results of simulations using this kernel are shown in figure [7] along with the predicted 
average high and low densities from the free energy. To calculate the coexisting densities 
of the high and low density sites we form a double tangent construction on the free 



energy 35 . As predicted the finite range nature of the interactions effectively creates 
a surface tension, but does not significantly alter the coexistence densities. Further, as 
expected we see no dependence in the steady state on our choice of Vq and a. 



5.2.2. Anisotropic Kernels For anisotropic coarse graining kernels, however, the 
situation is more complex. One simple and natural way to introduce an anisotropic 
kernel is to account for the finite volume of bacteria by stating that there can be at 
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Figure 8. Free energy density of the exclusion model where the occupation of a site 



is constrained to be smaller or equal to 100 particles by the choice of rates (81) 



most rim bacteria on each lattice site and taking the occupancy at the arrival site as our 



n in equation (43). This forms a generalisation of the partial exclusion process 25 ,36 



and results in the jump rates 

dt{nt,nT) = dt[l- — ]. (81) 



rir 

In this case the effective free energy is limited to a region where p < Pm and this 
section is sketched in figure |8} A double tangent construction amounts to looking for a 
density piow such that the tangent of the free energy density at this point meets with the 
free energy density at p = pm, as can be seen from inspection of figure |8] This amounts 
to finding p such that 

^-l)=log^, (82) 

Pm ^ Pm 

which can be solved numerically and yields for the low density piow that coexists with 

Pm 

— = .203188. (83) 

Pm 

For a total average density larger that piow, the stable phase should thus be a 
combination of two phases, one with p = piow and one with p = pm, the ratio of 
the amounts of the two phases being set by constraint on the global mass. 

Interestingly, although the theory correctly predicts a change from a flat profile 
to phase separation, on examining the results of simulations of the underlying lattice 
system we found that the densities into which the system separates do not correspond 
to those predicted by the continuum theory. Indeed, while the continuum theory had no 
dependence on the tumble rate a or the coefficient of the jump rate Vq, the simulations for 
an anisotropic kernel showed a strong dependence on the ratio of these two parameters. 
The lower and upper densities both varied considerably with the average run length 
r = vo/a, as shown in figure |9| and below r = 4 we see no separation at all. 

This discrepancy is not limited to the particular choice of anisotropic kernel we 
use as illustration above and is general to any anisotropic choice of K^{x). We have 
also conducted simulations with smooth anisotropic kernels without a hard limit on the 
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Figure 9. Snapshots of typical density profiles for an average run length of 100 sites 
(red), 20 sites (green) and 10 sites (blue). Data recorded sX t — 1000 with a = 1, 
n„j = 100 and from 10, 000 particles. 



number of particles per site, and saw exactly the same qualitative effect. Thus we have 
seen that the relevant factor is indeed the isotropy of the kernel but although the origin 
of the difference has been established, a comprehensive explanation for the variation 
between isotropic and anisotropic kernels is yet to be formulated. 

Note that it is not in itself surprising that a fluctuating hydrodynamics developed 
to describe smooth profiles fails to quantitatively amount for the coexisting densities 
of profiles alternating high and low densities. One of the reasons our fluctuating 
hydrodynamics work so well for the isotropic case is that we always consider large 
occupancies on each lattice site. This means that the model is close to mean-field in 
the same sense as the large spin limit of a spin chain is well described by a continuous 



spin chain 25 . Smaller occupancies would lead to quantitative differences between the 



predicted coexisting densities and those predicted from the fluctuating hydrodynamics. 



even for isotropic kernels. Furthermore, the Ito drift that was neglected in equation ( 75 ) 



for symmetry reasons would not vanish for the anisotropic case. In fact even for the off- 



lattice model, the fluctuating hydrodynamics developed previously 13 is only valid for 
isotropic kernel and the quantitative mismatch between the fluctuating hydrodynamics 
and the simulations on lattice for anisotropic kernels are thus not that surprising. We 
nevertheless now try to shed some light on their origin. 

5.2.3. Stability Analysis One way we can analyse the difference between the isotropic 
and anisotropic interaction kernels is to examine the dynamic stability of the two 
systems. We consider a one dimensional system evolving from a homogeneous state 
under a small perturbation and determine whether the system is dynamically stable or 
unstable, whether the perturbations will, on average, grow or shrink. 

One possibility is that the discrepancies we saw in figure |9] between the theory 
and simulations arise from the assumptions behind the diffusive limit taken in the field 
theory, see section|4]and Appendix A We therefore start from the continuum microscopic 



mean field equations for homogeneous and isotropic rates, i.e. after the continuum limit 
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has been taken but before the diffusive hmit, 



_P_ 

Pm 



a p 



P 



vV 



P 



1 



_P_ 

Pm 



+ 



2 2 

a p+ a p~ 



2 2 

We expand around a flat profile and Fourier transform. We take p^{x) 
(5^ exp(z q x) and arrive at 



(84) 

(85) 
Po/2 + 



51. 



v6jiq 1 



Po 

Pm 



+ viq 



Po 



'm ^ 

+ + (87) 



Po 

Pm 



viq 



Po 



We can re-write these two equations in matrix form as 



5n 



-viq[l 



3po 

2/3m 



Viq 



po 

2pm 



+ - 



viq[l — 



3po 

2 pm 



5n 



5: 



The eigenvalues of this matrix, which will then tell us whether the flat profile is stable 
or unstable to small perturbations, are 

„2 / . \ /o . \ \ 1/2 



a 



A±(g) = --± 



a 



I 2 2 

+ V q 



Po 

Pm 



2po 



- 1 



9) 



4 " V PmJ \P 

It is clear that one of these eigenvalues will always be negative while the other is negative 
for Po < Pm/2 and positive for po > Pm/2. Hence a homogeneously flat profile is stable 
when the average total density is less than half the maximum density and unstable above 
that, with no dependency on run length. This corresponds precisely with the stability 



predicted by the free energy derived in section |5.2[ That stability analysis was derived 
from a free energy which considered only the total density and was itself calculated only 
after assuming a diffusive scaling. That the diffusive scaling does not alter the criterion 
for instability justifies taking that limit and implies that the discrepancy between our 
lattice simulations and continuum free energy arises from another factor. 

We turn, then, to consider the dynamic stability of the lattice dynamics directly. 
Beginning with the mean field equations for the anisotropic partial exclusion process. 



n- 



dnj 



dn. 



+ ^ (90) 

« I (l-—j+^ 2- 

no + J2q ^"q exp(2 q k). After some algebra, 
detailed in Appendix B , we arrive at a condition for there to exist positive eigenvalues, 
i.e. for a flat profile to be unstable to small perturbations. Specifically we see instability 
whenever the run length r = d/a satisfies the following inequality 

1 



we expand around a flat profile, taking 



r > 



2 1 



no 



2 no 



1 



(92) 



Lattice Models of Nonequilihrium Bacterial Dynamics 



24 



100 


1 1 


1 1 


80 






60 


1 


II 1 


40 






20 









1 1 1 1 



0.0 0.2 0.4 0.6 0.8 1.0 

Fractional Density 



Figure 10. A flat profile is stable when in region I, but unstable in region II, for a 
system with exclusion and homogeneous, isotropic jump and tumble rates. The x-axis 
is the fractional density, i.e. n/rim. 



Graphing this we can see that a flat profile will be stable in region I of figure 10 and 
unstable in region II. We see that for short run lengths the range of densities in which 
the system will spinodally decompose is restricted and at run lengths below 4 sites there 
is no separation at all. Conversely, in the limit that the run length tends to infinity, 
i.e. where we effectively have two totally asymmetric partial exclusion processes on the 
same lattice, the system is unstable for any density between — = 0.5 and — = 1. This 

' rim rim 

instability is in accordance with that seen in our simulations. 

In our simulations we found that when we replaced the anisotropic kernel in the 
interaction terms with an isotropic one we recovered the density profiles predicted by 
the continuous theory. We can also analyse the effect of an isotropic density kernel on 
the dynamic stability. 

Consider now the microscopic mean field equations 



n 



n 



dnU (l - ^ E - {^-^Y. ^>^-) - + (93) 



3 3 

In general can take any values, we enforce only that they are both normalised, i.e. 
that 'Y^- = 1 . Relaxing this constraint would effectively re- normalise the maximum 
density. This more general interaction reduces to the simple exclusion case if we take 
Kf = 5-1-1 J and to the zero-range case if we take = 6oj. 

For an isotropic kernel, where Kf = Ki, when we expand around a fiat profile we 
find that there exists a q such that A+(g) is greater than if and only if rim > ""-o > 0.5 rim, 
see Appendix B for details, which matches the condition we derived from our continuum 
free energy. Thus for an isotropic interaction kernel we recover the continuum stability 
result while an anisotropic kernel will, in general, not produce the same result. It 
thus seems the error comes from the continuum limit itself, which is not valid for 
anisotropic kernels. While the stability analysis accounts qualitatively for the difference 
between isotropic and anisotropic kernels a theory which quantitatively accounts for the 
differences in coexistence densities at steady state remains to be constructed. 
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Figure 11. Snapshots of the two dimensional system using the partial asymmetric 
exclusion and parameters n,„ = 20, = 10 and a ~ \. Sites with a density greater 
than 0.6rim are marked in red, while those with a lower density are left blank. Left: 
allowing only nearest neighbour hops and no diagonal movement. Right: allowing for 
diagonal hops. Both simulations performed using 400, 000 particles and recorded at 
t = 2500. 



5.3. Finite Range Interactions in 2D 

Most cases of physical interest require a model in more than one dimension; it is 
therefore natural to extend our analysis to higher dimensions. This can be done on 
lattice relatively easily. On a square lattice in two dimensions we allow the particles to 
jump between nearest neighbours. We therefore consider 4 types of particle now instead 
of 2. We find that in two dimensions, the behaviour of the run-and-tumble crowding 
model is qualitatively the same as in one. The system separates into regions of high and 
low density, where those co-existent densities are given by the same free energy as in 
one dimension. The field theoretic approach developed in section |4] indeed generalizes 
straightforwardly to higher dimensions and yields the same fluctuating hydrodynamics. 
We find, however, that allowing only nearest neighbour hopping results in an unrealistic 



surface tension because of the anistropy of the lattice 37 , 38 ; the regions of high and 



low density form elongated, and thus anisotropic, domains, see figure 11 (left). 

To correct this unphysical characteristic of our model we extend our dynamics 
to allow next to nearest hopping along diagonal directions and consider 8 species of 
particles, where the jump rates in the diagonal directions are scaled by a factor of a/2. 



The droplets then coarsen into more realistic curved domains, see figure 11 (right). 



The stability conditions remain the same as do the coexistence densities and the 
discrepancy between isotropic and anisotropic kernels. We examined simulations with 
both the anisotropic partial exclusion kernel, and the isotropic, Gaussian kernel on 
lattices of 500x500 sites for densities above the spinodal point. The systems were seen 
to separate into droplets of higher and lower density which then coarsened into discrete. 



contiguous domains, see figure 11 



We measured the coarsening of these domains and found them to scale as t^/^, 
as we would expect for conserved model B type dynamics 39 . An approximate 
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jj^g Distance (re-scaled) 

Figure 12. Left: The average diameter of domains over time for an anisotropic 
kernel (green), Gaussian kernel with fc = 1 (red) and k = 2 (blue). The solid lines 
show asymptotic behaviour with an exponent of 1/3. Right: The correlation function 
C{j, t) for the anisotropic kernel simulations. The x-axis has been rescaled x — ?> x/t^-^^ 
so that data computed at t = 250 (red), t = 500 (blue), t = 1000 (green), t = 2500 
(black) and t = 5000 (magenta) superimpose. Both figures derived from data for 
240, 000 particles and with parameters Um = 10, a = 1 and do = 10 on a lattice of 
200 X 200 sites. 



measurement of the size of the domains was calculated by randomly sampling the system 
and measuring the horizontal and vertical size of the encountered droplet at that point. 
Mathematically, we define 

Lxihj) = max {A; G N : \nij — rii^mjl < n*, Vm G [0, k]} 

+ max {A; G N : \nij — ni_m,j\ < n*, Vm G [0, k]} , (95) 

where n* is an arbitrary cutoff to distinguish the two domains but ignore random 
fluctuations. Computations were run with a number of choices for n* and the particular 
choice of cutoff was found to have no significant effect on the results. We calculate the 
vertical size in an analogous fashion and average the lengths over a large number of 
points on the lattice. Though this does not give an exact measure of the droplet size 
it is sufficient to show the scaling of the domain size with time whilst being quick to 
calculate numerically. 

Using this procedure we determine that the domains increase in size with an 
exponent of approximately 1/3, i.e. (L^) (f ) = {t' /ty^^ {L^) (t). The two-point, 
connected, equal time correlation function C{j,t) = {ni{t) rii+jit)) — {riiit)) {rii^jit)) 
was also calculated numerically from the data and fit reasonably with a re-scaling 
C(x, t) = C{x/a^^^, t a). The data for both these measurements can be seen in figure 
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As can be seen from the figures the choice of kernel does not change the coarsening 
exponent, only the relative speed of coarsening, with interactions over a larger number 
of sites taking longer to reach a steady state than those with shorter ranges. 

6. Conclusion 



In this work we have presented several classes of lattice models based on the run-and- 
tumble dynamics of certain species of bacteria, notably Escherichia coli. We calculated 
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the exact steady state probability distributions for both inhomogeneous, anisotropic, 
non-interacting and zero-range interaction models. For more general types of interaction 
we used a field theoretic approach to derive the continuum fluctuating hydrodynamics 
and from there derived a mapping to a free energy like functional describing the steady 
state profile for a crowding interaction. We analysed the linear stability of both the 
continuum and lattice microscopic mean field equations and isolated a condition on the 
coarse graining we employed in our interaction terms for the zero-range free energy to 
work as a mean field theory. 

Our work builds on earlier treatments of run-and-tumble bacteria where interactions 



between bacteria were not included 12 , 40 . It provides a lattice counterpart to prior 



continuum approaches and qualitatively reproduces results obtained off lattice, where 



a similar, though not exactly equivalent, density dependence was considered 13 . Our 



approach on lattice provides a microscopic justification for manner in which this density 
dependence manifests; previously it was added in an ad-hoc manner after the diffusive 
approximation had been taken where in this work the dependence is intrinsic from the 
microscopic definition of the dynamics. That we produce qualitatively similar results 
justifies the way in which the dependence was inserted there. Our work also reveals 
a condition on the way in which the coarse graining in this density dependence must 
be performed for free energy mapping to work; it must be taken isotropically, so that 
particles moving in any direction feel the same effective local density. 

This work provides a means to simulate microscopic run-and-tumble dynamics 
efficiently, which is particularly important in two or more dimensions where microscopic 
simulations off lattice are very computationally expensive. It illustrates potential 
hazards in comparing lattice simulations and continuum theoretical predictions and 
offers some insight into how to avoid those problems by carefully choosing how to 
implement non-local interactions in the lattice dynamics. Note that early studies of 
off-lattice run-and-tumble dynamics actually had to discretise space to do simulation 
for position dependent swimming speed in order to compare with their theoretical 



predictions 40 . We hope the techniques we presented in this paper will develop the 
use of lattice simulations in the bacterial context much further, for example, to study 
pattern formation. 

We hope that the lattice approach will provide new exact results for run-and-tumble 
dynamics in higher dimensions, generalizing previous studies to more general external 
potentials for sedimentation and trapping 
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Last, there are some open question regarding the coexistence densities for 
interacting bacteria. Although our fluctuating hydrodynamics correctly predicts the 
existence of phase separation, the coexistence densities are accurately predicted only 
for isotropic kernels and large lattice occupancies. Quantitative predictions beyond this 
case are yet to be derived. More generally, explicit two body interactions and more 
general non-linear interaction kernels have not been looked at yet. These are a few of 
the outstanding questions among the many concerning this new and interesting class of 
models. 
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Appendix A. Hydro dynamic Limit and Scaling of Fields in the Action 

The microscopic action for the non-interacting, homogeneous and isotropic model can 
be written as 



-(pi+l—pi+Ji+l—Ji) _j_ ^Pi+l—pi — {Ji+l—Jt 



h) _ 2^ 



Ji 



_|_ _Z. [ ^-{Pi+I-Pi + Ji+I-Ji) _ Qpi+l-pi-{Ji+l-Ji) 



2 



) 



+ ^{Pi+i - pi- 



d 



+ 



(A-1) 



The continuous hmit can be taken by exphcitly introducing the lattice spacing a and 
making the substitutions 



Pi^ap{x); pi^p{x); d^va^; 

i 

Ji J{x); Ji J{x); Vi^ aV + Ja^A 



l=La 



dxa ^; 



(A-2) 



where Vj is the discrete gradient, e.g. Vipi — Pi+i — Pi. This overall substitution and 
the Taylor expansion of the action then gives 







dt / da; 







pp + v~^jj - vpVJ - JVp + ^ (e^-' + e 



2J 







iv 



(A-3) 



where the neglected term 5*1 is given by 

Si^ - [ dt [ dx \vp[ - JaJ + (Vp) V2 + (VJ) 72] + J[ - Ap/2 + VpVJ] 

+ v/2{Vp - VJ/v){Vp - VJ)j + 0{a). (A-4) 

To calculate the correct manner in which to rescale our fields let us begin by considering 
a system i times larger and rescaling t ti"", x ^ xi, p ^ p/i so that the action is 
given by 



dt I dx 





a 







pp + £—- i^'-^vpVJ - rJVp -e'^^pi + e-'-^ - 2 
V 4 



Av 



(A-5) 



For the J terms to not blow up we need to have J small. We therefore expand the 
exponentials to give 



S 



^- dt 

Jo Jo 



dx 



JJ 



a 



pp^t r-VVJ-^JVp-rapJ^ + r+'-JJ . (A-6) 

V V \ 
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If we explicitly take J to scale as J — )■ J£ 



S 



dt / dx 







PP + 



fll-l3-5'_ 



we get 



(A-7) 



and J Ji-^ 
■M 

V 

- r-^JVp - r-^^apj^ + r-p-^^ 

Now, we need the coefficient of each term to be of order 1 or smaller so that no terms 
blow up so 

1-/3-5 < 0; a-l-(3 < 0; a-S < 0; a-2/3 < 0; l+a-l3-S < 0.(A-8) 

However, as we do not want to simply be left with p = we need a — 6 = and as 
we also want to retain a noise, which corresponds to the term, we also require that 
a — 2/3 = 0. That leaves our action 



S 



dt / dx 



'^vpVJ - JVp - apJ^ 



(A-9) 



Which tells us that /3 < 1, /3 > 1/3 and /3 > 1 which imply that /3 
and 5 = 2. Injecting these scalings back into the action gives 



1 and hence a = 2 



S = — dt dx pp + 

Jo Jo 

where the macroscopic observation time r 



vpVJ - JVp + ap.P + 



aJJ 



(A-10) 



T/i"^ is supposed to be of order 1. The 
term in jj is thus irrelevant and we can check that the hydrodynamic action 



So 



dt / dx 



pp — vpVJ — JVp + apJ + 



aJJ 



(A-11) 



is invariant under further diffusive scaling. Note that the scaling of the fields considered 
here is arbitrary and we could choose to look at currents J, J larger than l/£^, This 
would correspond to trajectories whose probabilities are smaller than exp(— £), which 
are not correctly described by fiuctuating hydrodynamics and large deviations. One 
can also check that under this rescaling, the action Si stays of order 1 and aSi is thus, 
indeed, negligible. 

For one isolated bacterium the run-and-tumble dynamics is a variant of a random 
walk and is known to be diffusive at large scales [l2, 13 . It is therefore not surprising 
that we find a = 2. In the presence of interactions (as in section 4.2), a uniform density 
profile of bacteria will continue to exhibit diffusive behaviour; the interactions will simply 
rescale the diffusivity. If interactions cause the profile to become unstable, however, the 
model can nevertheless give rise to length scales which grow in a non-diffusive manner. 



as for example for coarsening (see section 5.3). 



Appendix B. Stability Analyses 



Beginning with the mean field equations for the partial exclusion-like interaction. 



Lattice Models of Nonequilibrium Bacterial Dynamics 



32 



we can expand around a flat profile and take rif. — uq + 5^ exp{iqk) to investigate 
the linear stability. In matrix form the resulting equations can be written as 



\ rim J ^ 
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where = {6^,6^ ) as before. Defining the run length as the ratio d/a = r, we can 
write the eigenvalues of this matrix as 



A±(g) = a 
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Again, one eigenvalue is always negative while one can be positive or negative. In 
this case, however the condition for stability is no longer independent of q and r. In 
particular, we find that A_|_ > when 
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no 



-2r^ + 4r^ 
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Further, if wc want to know only under what conditions on r and — there exist any 

Tim 

positive eigenvalues, we can consider the simpler condition 
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For the more general case where the microscopic mean field equations are given by 
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we can perform a similar analysis. Once again we expand around a fiat profile, this time 
to obtain 

^ a 
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Now, we can define a new function ^^{q) = -^j" 6xp(2 j g). Witli tliis definition, and 



the normalisation of K^, we can simplify equations (A-ll). 



5. 



'd 1 



no 



dnp 
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For isotropic kernels, where n^{q) = K,{q), it turns out we can write the eigenvalues in 
a relatively simple form and, even without knowledge of the specific shape of Kj, we 
can analyse the conditions under which there will exist positive eigenvalues. We start 
by writing the eigenvalues of the matrix in equation (A-12) as 
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A±(g)=a( -2+nl 
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The larger of these two eigenvalues will be positive if 
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Further, as K{x) is always positive, K,{q) will have a maximum at g = 0, where k(0) = 1, 
and so we can examine the simpler condition 



2r' 



r + 



2rnQ Ar'^nl Gt'^Uq 



> 0. 



(A-15) 



When this inequality is fulfilled we will see instability. Given that r must be positive 
this means any isotropic density kernel will be unstable in the range uq G [0.5 rim, nm]- 



